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Abstract 

The paper addresses a numerical method for solving second order elliptic partial differential 
equations that describe fields inside heterogeneous media. The scope is general and treats the 
case of rough coefficients, i.e. coefficients with values in L°°(fi). This class of coefficients 
includes as examples media with micro-structure as well as media with multiple non-separated 
length scales. The approach taken here is based on the the generalized finite element method 
(GFEM) introduced in [5], and elaborated in [3], g] and [25]. The GFEM is constructed 
by partitioning the computational domain Q into a collection of preselected subsets = 
1,2, ..m and constructing finite dimensional approximation spaces over each subset using 
local information. The notion of the Kolmogorov n-width is used to identify the optimal local 
approximation spaces. These spaces deliver local approximations with errors that decay almost 
exponentially with the degrees of freedom Ni in the energy norm over Wj. The local spaces 

are used within the GFEM scheme to produce a finite dimensional subspace S N of H 1 (Q) 
which is then employed in the Galerkin method. It is shown that the error in the Galerkin 
approximation decays in the energy norm almost exponentially ( i.e., super-algebraicly) with 
respect to the degrees of freedom N. When length scales "separate" and the microstructure is 
sufficiently fine with respect to the length scale of the domain cji it is shown that homogenization 
theory can be used to construct local approximation spaces with exponentially decreasing error 
in the pre-asymtotic regime. 
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1 Introduction 

Large multi-scale systems such as airplane wings and wind turbine blades are built from fiber rein- 
forced composites and exhibit a cascade of substructure spread across several length scales. These 
and other large composite structures are now seeing extensive use in transportation, energy, and in- 
frastructure. The importance of accurate numerical simulation is ever increasing due to the high cost 
of experimental testing of large structures made from heterogeneous materials. The computational 
modeling of such heterogeneous structures is a very large problem that requires the use of parallel 
computers. In order for a numerical method to be adequate it must be able to utilize many local 
computations performed independently on single processors or clusters of processors of reasonable 
size. The approach taken here is based on the Generalized Finite Element Method (GFEM) intro- 
duced in [5] , and elaborated in [3] , [4] , |25j . It is a partition of unity method [5] which utilizes the 
results of many independent and local computations carried out across the computational domain. 
The GFEM is constructed by partitioning the computational domain fi into to a collection of prese- 
lected subsets uji, i = 1,2, ..m and constructing finite dimensional approximation spaces ^ over each 
subset using local information. The specific way in which the partition is carried out is special to this 
method [5] , [3j and the details are discussed in section [2j Since each space "I^ is computed indepen- 
dently the full "global" solution is obtained by solving a global (macro) system which is an order of 
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magnitude smaller than the system corresponding to a direct application the finite element method 
to the full structure. The GFEM approach provides an opportunity for the significant reduction of 
the computational work involved in the numerical modeling of large heterogeneous problems. 

In this article we show how to achieve optimal accuracy within the GFEM approach. The key 
point is to note that the approximation error of the GFEM is controlled by the corresponding 
approximation error of the local approximation spaces i = l,...,m, see section [2j Therefore 
the goal is to identify optimal local approximation spaces. Our approach is naturally guided by the 
notion of the Kolmogorov n-width [22* which measures the ability of an increasing sequence of finite 
dimensional subspaces of a prescribed Banach space B to approximate any element inside B, see 
section [3] Using the solution of the spectral problem associated with the n- width we are able to 
identify a new class of approximation spaces We show that these finite dimensional spaces are 
able to approximate the solution on u>i with errors that decay almost exponentially with the degrees 
of freedom Ni in the energy norm. The overall method for constructing the local approximations is 
general and applies to subdomains belonging to the interior of the computational domain as well 
as those that intersect the boundary of the computational domain. The optimal local approximation 
spaces are identified for interior sub doma ins in section [3~T] and for those touching the boundary of 
the computational domain in section [3. 2 1 

The optimal local spaces are then combined within the GFEM scheme to produce a finite 
dimensional GFEM subspace S N of i/ 1 (SI) which is employed in the Galerkin method. It is shown 
that the corresponding error in the Galerkin approximation decays in the energy norm al most ex- 



ponentially ( i.e., super-algebraicly) with respect to the degrees of freedom N see section 3.3 In 
section [4] we discuss the main issues involved in the implementation as well as estimates of the 
computational work associated with this numerical approach. 

We show how to construct nearly optimal local approximation spaces using homogenized coeffi- 
cients when the subdomain u>i is sufficiently large with respect to the length scale of the heterogeneity. 
The homogenization limit for the n-width and the optimal approximation space is identified in sec- 
tion [5] This identification is established within the general homogenization context described by 
i/-convergence and G-convergence, [T7] , [53] . These results are applied to heterogeneous media with 
micro-structure that has uniformly fine variation with respect to the length scale of the domains 
LJi. Here a uniformly fine microstructure is defined to be one that can be identified as belonging 
to a sequence of microstructures characterized by a sequence of length scales e = \, k = 1,2,... 
and coefficients {A e } €> o that converge to a homogenization limit described by a matrix of constant 
coefficients A. For this case we provide examples that illustrate how to construct local approxima- 
tion spaces with errors that decay exponentially in the pre-asymptotic regime, see section [6j The 
examples corroborate the exponentially decreasing error observed in the numerical simulations for 
finely mixed dispersed inclusions carried out in [25] . 

The homogenization theory developed here provides motivation for some rules of thumb for 
choosing the size of subdomains tjj in the implementation. Here the size is chosen large relative to 
the local length scale of the heterogeneity but small enough such that the heterogeneity is statistically 
uniform within it. The specific details are presented in section [6] 

We conclude noting that there is now a large and rapidly growing literature devoted to the 
numerical analysis of multi-scale media. Several contemporary mathematically based approaches 
include the Multiscale FEM [TJ], [13], global changes of coordinates for upscaling porous media 
flows [TO], [H] , the heterogeneous multiscale methods (HMM) [EMS], [H], an adaptive coarse scale 
- fine scale projection method [18], numerical homogenization methods for L°° coefficients based 
on harmonic coordinates and elliptic inequalities [19], [20], [6J, subgrid upscaling methods [1] and 
global Galerkin projection schemes for problems with L°° coefficients and homogeneous Dirichlet 
boundary data |15| . 

1.1 Problem formulation 

Let 11 e l d be a bounded domain with C 1 smooth boundary dfl. In this article we consider the 
elliptic differential equation 

- div(A{x)Vu(x)) = f(x), VieO (1.1) 
with either Neumann boundary conditions prescribed on the boundary <9S! given by 

n ■ AVu{x) = 5 , xGdfl (1.2) 
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where n is the outer unit normal vector or Dirichlet boundary conditions 



u(x) = q, x 6 (9f2. (1.3) 

In forthcoming work we will address the case of non-smooth boundaries and the case where both 
Dirichlet and Neumann boundary conditions are prescribed on different parts of the boundary. 

We assume that A(x) is a d x d symmetric matrix with measurable coefficients ciij(x) g L°°(f2) 
satisfying the standard coercivity condition 

a*(x) | v | 2 < v*A(a;)v < 0*(x)\ v | 2 ,Vx G Q (1.4) 

where v G R d is an arbitrary vector and 

< a < a*(x) < /3*(z) < /3 < oo, Va: € O. (1.5) 

For future reference we will denote this class of coefficients by €. Here we will suppose that / G 
H k (Q), k > 0, g £ H~ 1/>2 (dQ), together with the consistency condition J g(1 gds + J Q fdx = for the 



Neumann boundary condition (1.2) and that the Dirichlet data q is an element of H^ddQ). The 
weak solution belonging to u G H (O) exists and for Dirichlet boundary conditions it is unique and 
for Neumann boundary conditions is unique up to additive constant. 

In this work we investigate the accuracy of the GFEM for the approximation of the exact solution 
uq. Here the objective is to find an approximate solution u T G H (Cl) for which 

II uq - u T ||f(n)< t 

where 

II u \\s{n)= (J AVu ■ Wudx 

is the energy norm. In this treatment we address the scalar problem noting that the ideas used in the 
approach presented here apply with out any modification to second order elliptic systems including 
the system of linear elasticity. 



1/2 

(1.6) 



1.2 A typical example: fiber reinforced composites 

For the purposes of this article we have chosen to work with coefficient matrices belonging to L°°(f2) 
subject to standard coercivity and boundedness conditions. This choice reflects our intention to 
describe generic situations for which there can be several non-separated length scales of variation 
inside the heterogeneous media. In this treatment we shall assume that the coefficient matrix 
describing the media is known. To fix ideas we discuss the problem of determining the coefficient 
matrix A associated with a sample of fiber reinforced composite material described in [2J. A fiber 
reinforced material consists of two components the fiber and the host material commonly referred 
to as the matrix. The material shown here is taken from the center of a composite plate of 36 plies, 
divided into 9 groups each containing 4 plies. The orientation of the fibers alternates between 0° and 
90° from group to group. The plys are HTA/8376 unidirectional prepreg fiber composites produced 
by Ciba-Geigy. The sample is a rectangular plate of length 300mm and width 140 mm. The nominal 
ply thickness is 130 fim. Figure [T] shows the cross section of a group of four plies consisting of 16275 
fibers. The cross section of each fiber is roughly circular with a fiber diameter of about 7/xm. It is 
clear from Figure [T] that the material coefficients have variation across several length scales. These 
include "matrix rich" zones between the plys as well as variation in fiber alignment between groups 
of plys. At the smallest length scale Figure [2] shows that the coefficients are piece wise constant 
taking one value in the fiber and a different value in the matrix. This composite sample has been 
mapped in the study [SJ and the relative fiber positions are known exactly and so the coefficient 
matrix can be determined. In general it is impossible to record the relative location of every fiber 
for an entire structure made from a composite material. Hence some stochastic information needs 
to be extracted from the structure and used to describe the coefficient matrix. An investigation of 
the different stochastic data that can be obtained from composite samples is taken up in [2]. As 
mentioned earlier we will assume that the problem is deterministic with well defined coefficients. 
Future work will address the problem of determination of the approximation error for stochastically 
defined coefficients. 
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Figure 1: Fiber-reinforced composite. 




Figure 2: Microstructure. 



2 The generalized finite element method (GFEM) 

The GFEM is based on the partition of unity method (PUM) and is introduced in [5] as a method 
for the numerical solution of elliptic PDE with rough coefficients. The method is further elaborated 
and extended to other application areas in the works of [3] , [4] , [24] and [25] . We briefly summarize 
the main ideas and results of GFEM. For more details see [3]. We recall the scalar Dirichlet and 
Neumann boundary value problems defined in section For the Dirichlet problem define the 
hyperplane H%(Sl) = {it G | u = q{x) on <90} and set H&(fl) = ff*(0) for q = 0. The weak 

solution of the Dirichlet problem uq £ H q (ft) satisfies 

B(u ,v) = F(v), (2.1) 

for all v £ Hq(£1) where 

B(u,v) = / A(x)\7u ■ \7vdx, and F(v) = / fvdx. (2.2) 
Jn Jn 



For the Neumann boundary value problem the solution uo € ii^(O) satisfies (2.1) for all v € 
£^(0) and 



F(v) = / /vcfa; + / gvds 



in Jon 

For future reference we note that the energy norm (1.6 1 is given by || u ||,£(m= {B(u, u)) 1 / 2 . 



We now recall that for a N dimensional space S(ft) C H 1 (O) that the associated Galerkin solution 
of the Neumann problem u s G S(n) C iif 1 (f2) satisfies B(u s , v) — F(v) for all v £ 5(0). For a fixed 
tolerance e > if there is a ip € 5(0) such that || Uq — tp ||g(n)< £ then it is clear that the Galerkin 
solution satisfies |j uq — u s || £ (si)< e. A similar conclusion holds for the Galerkin solution of the 
Dirichlet problem. Let So(0) be an N dimensional subspace of i?o(0) and set S q (il) = Sq(CI) © $o 
where $0 is a particular function belonging to iJ^(O). Then the Galerkin solution of the Dirichlet 
problem u s £ S q (fl) satisfies B(u s , v) = F(v) for \/v £ Sq(£1). Again it is clear that if there exists a 
tp £ Sq(Sl) such that || uo — i/j ||£ ( n)< £ then || uq ~ u s \\ £ ^< e. Hence the main task in constructing 
a Galerkin numerical solution is the selection of S(Cl) for the Neumann problem and the selection 
of the hyperplane S q (Sl) for the Dirichlet problem. 

We introduce the Galerkin approximation delivered by GFEM and discuss the associated approx- 
imation error. Since most physical situations are described by Neumann boundary conditions we 
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address this case and note that the Dirichlet case follows identical lines. Let {Oi} 7 ^ 1 , be a collection 
of open sets that cover the computational domain, i.e., ft C U^Oj, and we introduce the partition 
of unity subordinate to the open cover denoted by 4>i € C 1 (Oi), i — 1, ...m. We relabel interior 
sets Oi C ft as u>i — Oi and for sets Oi that intersect the boundary of ft we write Ui = Oi ft. 
Here U^LjWj = ft and we assume that each point x £ ft belongs to at most k subdomains w;. The 
functions <f>i, i = 1, . . . , m have the following properties 



< fa < 1 

</n(x) = 

m 

E<m*) = 1 



max | fa{x) 
max | Vfa(x) 



i = 1, ... ,771, 

for x € i — 1, ..to 



< Ci, i = 1, ...m 



< 



diam(tUi) ' 



z = l, ...m, 



(2.3) 
(2.4) 

(2.5) 

(2.6) 

(2.7) 



where the constants C\ and Ci are positive and bounded. Here diam{ui) denotes the diameter of 
uii. 

Next we introduce local approximation spaces associated with each tOj. For this case let ^i be a 
finite dimensional subspace of i? 1 (wi) of dimension ATj. The trial and test spaces for the GFEM are 
constructed from the local approximation spaces and are defined by 



Ni 



S(Si) = { V; ff> = E E ^ where e * 



(2.8) 



i=i j=i 



Note that although ^ only belongs to i/ 1 ^) the PUM construction ensures that S(ft) C i? 1 (0). 
Here the dimension iV of S(f2) is given by N = YliLi Ni- 

An analogous approach using local approximation spaces is used for the Dirichlet problem. Here 
the only difference is when UjOdft ^ 0. For this case if q = then fyj c i/ 1 (w J )niJQ (ft) and if q ^ 
then \tj C {^(iOj) n ifo(^)) © $j where $j is a particular function belonging to H l (u}j) n i?i(f2). 

It is now clear from the formulation that the approximation error of the Galerkin numerical 
solution for the GFEM is tied to the accuracy of the local approximation spaces. With this in mind 
we state the following approximation theorem [4] for the Neumann problem. 



Theorem 2.1. Suppose that there exists Q € i = 1, 



, m with 



I «o - Ci Hljo*) = (/ /8*(uo-Ci) 2 <fa0 1/2 <ei(0 

II wo - Ci IU(a,») < 



(2.9) 
(2.10) 



and Zet 



C(x) = ^2&(x)fa(x) 
»=i 

tfien f(ar) G iJ 1 (fi) and 

m 

||«o-C||xs ( n) < ^(E^iW) 2 ) 172 

i=l 

IK-CIM < (c 2 2 E(^«M^ 2 (^)) + Ci 2 E £ 2«) 1/2 ; 



(2.11) 
(2.12) 



where C\ and C2 are as in (2.6) and (2.7) respectively. 

Moreover if all local spaces U*,; contain the subspace of con stant functions then, for an appropriate 
choice of constant, the first term on the right hand side of (2.12) can be omitted as it is majorized 
by the second term. 
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We remark that Theorem |2.1| also holds true for Dirichlet boundary conditions. For this case 
the only modifications are the ones previously discussed for the subdomains ojj that intersect the 
boundary. Theorem |2.1| shows that the proper selection of the local approximation spaces are 
essential for obtaining optimal accuracy. In the next section we identify local spaces ^ that deliver 
a nearly exponential rate of convergence with respect to the dimension of the trial space. 



3 Optimal local approximation spaces and nearly exponen- 
tial upper bounds on their accuracy 

We identify the optimal local approximation spaces for GFEM and provide an upper bound on the 
accuracy of their approximation. We w ill consider the case when the exact solution uq solves the 



Neumann problem with / = in (1.1). During the course of the exposition we will indicate the 
modifications needed to treat the corresponding Dirichlet problem. The generalization to the case 
/ 7^ will be addressed as part of the implementation given in section [4] 

We recall the local approximation spaces introduced in section 2 denoted by \£j defined on the 
sets oji. To fix ideas we will assume that oji are the cubes of a given side length surrounded by a 
larger cube oj*. We will distinguish two cases depending on if the set Wj, lies within the interior 
of O or if uji n dfl ^ 0. It will be shown that the overall approach to constructing optimal local 
approximation spaces for these two cases is the same. We drop subscripts and consider concentric 
cubes oj C oj* with side lengths given by a and a* — (1 + p)a respectively. In order to introduce 
the ideas we suppose first that oj lies in the interior of Q so that oj C oj* C £l. The energy inner 
products associated with these subsets are defined by 



[u, u)f ( w ) = / AVu-Vvdx ( u , v )£(u*) = / AVu-Vvdx. (3-1) 

J OJ J D* 

We shall utilize oj* to construct a finite dimensional approximation space over oj. For any open 
subset S of the computational domain f2 we introduce the space of functions Ha{S) defined to be 
the functions in H 1 ^) that are A-harmonic on S, i.e., v £ H 1 (S) and 

B(v,<p) = 0, y^eC^(S). 

Here Ha{oj) and Ha(oj*) contain local information on the heterogeneities and will be used in the 
construction of the optimal local basis. We introduce the quotient of Ha(oj*) with respect to the 
constant functions denoted by Ha(oj*)/R. It is clear that the solution u lies in this local space 
modulo a constant. 

In this method we choose to approximate elements in the s pace o f functions Ha(oj*)/M. restricted 
to oj. This choice is motivated by the Caccioppoli inequality ( |3.11[ ), proved in Appendex A, which 
is used to estimate the energy norm over oj in terms of the L 2 norm over oj*. Let P : Ha(oj*)/H — > 
Ha(oj)/M. be the restriction operator such that P(u)(x) — u{x) for all x € oj and u € Ha(oj*)/M.. 
The operator P is compact, this follows from the Caccioppoli inequality, and the compactness proof 
is given in the Appendix, see Theorem |A.1| 

Now we approximate by "n" dimensional subspaces S(n) C H\(oj)/M.. The accuracy of a par- 
ticular increasing sequence {S(n)} n=1 of local approximation spaces is measured by 

d(S(n),oj)= sup inf l|P "7 Xll£M - ( 3 - 2 ) 

u£H A (uj*)/RX£S(n) \\U\\ £ ^,) 

A sequence of approximation spaces S(n) is said to be optimal if it has an accuracy d(S(n), oj) that 
satisfies d(S(n),oj) < d(S(n),oj), n — 1,2,..., when compared to any other sequence of approxi- 
mation spaces S(n). The problem of finding the family of optimal local approximation spaces is 
formulated as follows. Let 

A i , • f ■ r \\ Pu -X\\£(u) , - 

d n (oj,oj ) = mi sup ml — r — r . (3-oj 



G 



Then the optimal family of approximation spaces satisfy 

d n (w,o;)= sup inf — —r . (3.4) 

u£ffA(u')/EXe*»(w) 

The quantity d n {uj,u*) is known as the Kolomogorov n-width of the compact operator P see, [22j . 

The optimal local approximation space \E r rj , (cu) for GFEM follows from general considerations. 
We introduce the adjoint operator P* : H a (uj)/M. — > H A (u>*)/M. and the operator P*P is a compact, 
self adjoint, non-negative operator mapping Ha(uj*)/M. into itself. We denote the eigenfunctions and 
eigenvalues of the problem 

P*Pu = Xu (3.5) 

by {<fi} and {A^} and apply Theorem 2.2, Chapter 4 of [22] to see that the optimal subspace is 
given by the following theorem. 

Theorem 3.1. The optimal approximation space is given by ^ n (uj) = span{?pi, . . . , ip n }, where 
ipi = Ptfi and d n (uj,uj*) = ^X n +i- 

For the case considered here the definitions of P and P* show that the optimal subspace and 
eigenvalues are given by the following explicit eigenvalue problem. 

Theorem 3.2. The optimal approximation space is given by ^/ n (uj) = span{ip\, . . . , ip n } where 
ipi = Ptfi and tpi and Xi are the first n eigenfunctions and eigenvalues that satisfy 

(<Pi,S) eM = Aifo>i,<y) e( o> WeH A (u*)/R. (3.6) 



Proof. The eigenvalue problem (3.5) is given by 



(P*Pu,S) £ ^ } = X(u,8) £(u * ) , VS e H a (uj*)/R. (3.7) 



The theorem follows from (3.7| noting that 



(P*Pu, <5) £(0 = (Pu, PS) eM = (u, S) ew . (3.8) 

The next theorem provides an upper bound on the rate of convergence for the optimal local 
approximation. 

Theorem 3.3. Exponential convergence for interior approximations. 
For e > there is an N e > such that for all n > N e 

d n (u;,uj*) < e ~ n ^~'\ (3.9) 



The index N e is constructed explicitly in the proof of Theorem |3.3| given in the next section. Theorem 
3.3| shows that the asymptotic convergence rate associated with the optimal approximatio n sp ace 



we 



is nearly exponential for the general class of L°°(0) coefficients belonging to £. In section 3.2 
identify the optimal local approximation space for the case when u> touches the boundary of the 
computational domain. In that section we show that the convergence rate is also nearly exponential 
with respect to the degrees of freedom of the optimal local approximation space. We collect our 
results in section |3.3| and recall Theorem |2.1| to obtain the nearly exponential convergence rate for 
GFEM stated in Theorem El 
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3.1 Local approximation on the interior 



In this section we establish Theorem |3.3| To do this we construct a family of approximation spaces 
that exhibit nearly exponential convergence in the accuracy of approximation. The convergence rate 
for this family delivers an upper bound on the convergence rate for the optimal local approximation 
space described by Theorem |3. 2 1 

For future reference we introduce the decomposition of H 1 ^*) given by 

H^uj*) = H° a (uj*) + H^(uj*) +R. 

Here H^(uj*) is the subspace of Ha(uj*) given by the functions with zero average over uj*. The 
spaces H\(uj*) and Hq(uj*) are orthogonal with respect to the energy inner product (•, ^siu*)- The 
orthogonal projection from H 1 ^*) onto H^(uj*) is denoted by V A . 

The construction of the local approximation space is done iteratively. We start by introducing 
the the first n non-constant eigenfunctions (>,• 6 H 1 (uj*) of the Neumann eigenvalue problem 

(«i,iy)f( w ») = A,; / Viwdx, V w G H 1 ^*) 

J UJ* 

posed over oj*, i = 1, . . . ,n. The subspace spanned by these functions is denoted by S n (w*)- Next 
we introduce the span of A harmonic functions given by 

W n (uj*) = sp&n{wi £ H^(lu*) : Wi = Vi, on dui*, i = 1, . . . n} (3.10) 

One readily checks that W n (uj*) = V A S n (uj*). 

We define the family of approximation spaces !F n {ui,ui*) given by the restriction of the elements 
of W n {uj*) to uj. In what follows we first show that J-"„(oj, cj*) is a family of local approximation 
spaces with a rate of convergence on the order of n _1 / d , for d = 2,3. To show this we introduce 
a suitable version of the Cacciappoli inequality that bounds functions in the energy norm over any 
measurable subset O C uj* for which dist(dO, duj*) > S > in terms of the L 2 norm over uj* . 

Lemma 3.1. Let u be A-harmonic in uj* and belong to L 2 (uj*) n Hf oc (u:*). Then 

II u \\s (0 )< m 1/2 /6) II « h>w ■ (3.11) 



where (3 is defined in (1.5). 



The proof of Lemma |3.1| is given in the Appendix. Next we introduce the approximation theorem 
associated with the space W n (uj*) given by 

Lemma 3.2. Let u G H^(uj*) then there exists a v u G W n (uj*) such that 

l/d 

||«-«u|U»( w -)= jnf Alu-vh^KC^^a-^WuWe^,) (3.12) 

v£W n (uj*) V47T 



where a* is the side length of the cube uj* , jd is the volume of the unit ball in R d and C n = 
+ o(l)), ford =2,3. 

Proof. The lemma follows immediately from an upper bound on the quotient 

R = sup ml — 77— r s — -. (3.13) 

«eH»( u «)»ew»( u *) IrI|s(u*) 

Fix u G Ha(oj*)° and denote the projection of u onto W n {uj*) with respect to the energy norm 
II ' by 'P £ u. Choosing w — V £ u and noting that || (I — 'P £ )u\\ £ ( u] *) < IM|e(ui*) gives the upper 

bound 

R < sup -rr— r — - — -. (3-14) 

u£H%(uj')± W n (u*) \\ u \\s(uj*) 
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Since W n (u*) = V A S n {uj*) it follows that 

{u G H\{lo*) L V A S n (uj*)} = {ue H%(u*) J_ S n (u*)}, 
{u G £$(w*) -L S^w*)} C{«£ fl" 1 ^*) _L (5„(w*) + K)}, (3.15) 
where the _L in the second line of ( 3.15| ) is with respect to the L 2 (cj*) inner product. Hence 

D ^ IMU 2 (o;*) , NU a (u>*) 1 ,„ 

it < sup -tt— j| — - — i < sup -— j] — - — - = — — ■ (3.16) 

where is the largest Neumann eigenvalue associated with S' n+1 (w*). One has the elementary 

lower bound p n +i > Q^n+i where f n +i = 47rC~ 2 (<7*7y ,i )~ 2 is the corresponding Neumann eigen- 
value for the Laplacian with on squares (d = 2) or cubes (d = 3). The required upper bound on R 
now follows and the theorem is proved. 

Now we apply Theorem 3.1 to u — v u on uj C ui* and combine it with Theorem |3.2| to obtain the 
following convergence rate associated with the family of approximation spaces IF n {u),<jj*) given by 

Theorem 3.4. Let u G H\(uj*) then there exists an approximation v u G J- n (uj,UJ*) for which 



II u-v u \\ £ i u )= inf || u-v ||g (w ) < I(u),ui*)C n || u \\e( u *) (3-17) 



where 



1/d 1 + 

J(w,w*) = ^^(/3/«) 1/2 a«d C„ = n" 1 ^ + o(l)), d = 2, 3, (3.18) 

where is the volume of the unit ball in dimension d. 

Next we proceed iteratively to construct a family of local approximation spaces with a rate of 
convergence that is nearly exponential. For any pair of two concentric cubes Q C Q* we define 
FniQiQ*) to be the space given by the restriction of W n (Q*) on Q. We suppose that uj* is of 
side length a*. Let N > 1 be an integer and we suppose that oj is of side length a and a* = 
er(l + p). Choose uij, j = 1,2, ...N to be the nested family of concentric cubes with side length 
f (1 + p(N + 1 — j)/N) for which u — ujn+i C lon C ujv-i C ••• C wi = w*. We introduce the local 
spaces, J- n (uj,ujN), J- n (uj,ujN^i),...,J- n (uj,uji). Put m — N x n and we define the approximation 
space given by 

T(m,cj,LJ*) = F n {u,u}\) H h J"„(u;,cjjv). (3.19) 

The convergence rate associated with the local approximation space T{m, ui,U!*) is given in the 
following theorem. 

Theorem 3.5. Let u G if^(w*) and TV > 1 be an integer. Then there exists z u G T(rn, co, oj*) such 
that 

JV 



IU(a,)<? JV ||«IU( U .) (3-20) 



1/d 



and ? = 2 ^ l±a N{P/af' 2 C n 



' we 



Proof. In what follows we make the identification cj = and u;jv + i = w. From Theorem 3.4 
have that there exists v\ G J- n {bJ\,uj*) such that 

Vx\\s{ Ul ) <s\\u \\s(u") (3-21) 

Suppose next that for m = 1, . . . ,j there are functions v m G J>i(w m , CJ m _i) such that 

II ||f (w ,)<^ IMk"*) ( 3 - 22 ) 

m — 1 



9 



Applying Theorem 3.4 we see that there exists a Vj+i G J-" n (u>j+i,Uj) for which 



~ v m ) - v J+l \\ £[u]j+l) < q || u - V r, 



(3.23) 



and the induction step goes through. Choosing z u = ^ j v m delivers 



\\u- Z u \\e(u)< S N II U \\s(ui*) 



(3.24) 



and the theorem follows noting that z u belongs to 7~(m, oj, ui*). 

Next we make a choice for N. We choose N to be the largest integer less than or equal to ri 1 

for < 7. Thus m < n 1+ and to^ 1 < n and it follows that n a < m , and AT < mT +1 . On 

applying these inequalities we obtain 



q N < exp i —m~t+ 1 [ — In K 



lA*-7 
7 + 1 



In rn 



(3.25) 



1/d 



where K = 2 -^^(£) 1 ^ 2 (^^)- It is evident that decay occurs for the choice < 7 < | and 



(3.26) 



for m > N = (KeY 1+lS> ^ 1 ^ d J \ We set £ = dimcnsion{T(m, w, uj*)} and Theorem 3.5 together 
with (3.26 1 imply 



dg{ui,uj*) < sup 



inf 



for £ > N and Theorem |3.3| is proved. 



\u-x\\e(u) < e -£^b 



(3.27) 



3.2 Local approximation at the boundary 

Consider two concentric cubes C C uj* of side lengths a and a* = (1 + p)a respectively. We suppose 
that uj* n£7^ and w*fl30 ^ 0. The truncated cube u> is defined to be uj = Cnfi, and <9wn<9f2 7^ 0, 
see Figure[3j For this case we will assume that dfl is C 1 , i.e., the boundary can be represented locally 
as the graph of a C 1 function. The method presented here applies to both Dirichlet and Neumann 
boundary value problems. We will illustrate the ideas for the Neumann problem and make references 
to the Dirichlet problem when appropriate. Given a function u € Ha(Q) the goal is to provide a local 
approximation to u in uj. To this end we form a local particular solution u p given by the A- harmonic 
function that satisfies n ■ AVu p = g on uj* n dfl and u p — on dui* n ft. Writing u — u p + uq we see 
that n ■ AWuo = on uj* n dfl and uq — u on duj* n fl. If instead we have Dirichlet data then the 
particular solution u p satisfies u p = q = u on a;* fl dfl and u p = on duj* fl f2. The objective of this 
section is to find the optimal family of local approximation spaces that give the best approximation 
to i*o — u — u p in the energy norm over the set uj. To this end we introduce the the space of functions 
Hj\fi{fl n uj*) given by all functions v in ff 1 (17 n uj*) that are A-harmonic on fl n uj* and for which 

d v v = n ■ AVv = 

on dfl H u* . The analogous space of functions defined on uj is denoted by Ha.o(uj). Since we 
approximate functions with respect to the energy norm we introduce the quotient space of i?^,o(^n 
uj*) with respect to the constant functions denoted by HA,o(fl fl w*)/R. 

Now we introduce P : Ha.q {uj* D Q)/K — > Ha,o(uj) /M. given by the restriction operator defined 
by P(u)(x) — u(x ) for all x £ uj and u € H a ,o{uj* n 0)/R. The operator P is compact, this follows 



from Lemma 



A.2 



given in the Appendix. Let S(n) be any finite dimensional subspace of Ha,o{uS) /! 
and the problem of finding the family of optimal local approximation spaces is formulated in terms 
of the n-width of P. Let 

4(w,w*nfi) = inf sup inf ^f"~ X ^ M . (3.28) 

S(n) ueH A: „(ujT]n)/Rx£S(n) ||ti||g(u,*nn) 
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Figure 3: Truncated cube uj = C fl O. 



Then the optimal family of boundary approximation spaces {^ n (uj)} < ^' =1 for GFEM satisfy 

d n {cj,uj* nfi) = sup inf ll f"~ Xllg( " ) . (3.29) 

ueJ?A,o(w*nn)/Rxe*n(w) IMl£( w *nfi) 

Proceeding as before we introduce the adjoint operator P* : Ha(uj)/M. — ► Ha{uj* fl f2)/K and the 
operator P*P is a compact operator mapping Ha.o(uj* H fi)/R into itself. Similar arguments show 
that the optimal approximating spaces are given by the following theorem. 

Theorem 3.6. The optimal approximation space is given by ^>„(uj) = spanlipi, . . . ,i/j n } where 
ipi = Pifi and (p i £ i? J 4,o(< j -'*nJ7)/K and Xi are the first n eigenfunctions and eigenvalues that satisfy 

(<Pi,S) eM = x l (i Pt ,s) £( ^ nn) ,yseH Afi (io*nn)/R. (3.30) 



The next theorem provides an upper bound on the rate of convergence for the optimal local 
boundary approximation. 

Theorem 3.7. Exponential convergence at the boundary. 

For e > there is an N e > such that for all n > N t 

d^u^* nsi) < e -« w+1 ; . (3.3i) 



Theorem |3.7| shows that the asymptotic convergence rate associated with the optimal boundary 
approximation space is also nearly exponential for the general class of L°°(uj*) coefficients <£. 

We now give the prof of Theorem |3.7| The subspace of A-harmonic functions defined over fin a/ 
with zero mean is denoted by H A (il Dui*) and the subspace of elements belonging to Ha,o(Q Duj*) 
with zero mean over SI Do/ is denoted by H A (Slnw*). We introduce the L 2 (SI flu;*) norm closure 

of H A (fi fl uj*) denoted by H° A (flnw*). Useful properties of functions in H A (£1 n uj*) are listed 

below in Theorem 



3.9 



at the end of this section. It is shown there that v £ H A (S1 n uj*) implies 
that d v v = on d(l f 1 uj* and that v is A-harmonic in SI n uj*. 

Next we introduce the the first n non-constant eigenfunctions of the Neumann eigenvalue problem 
div(A VWj) = — XiVi over SI n uj* , i = 1, . . . , n. The subspace spanned by these functions is denoted 
by S n (Q n uj*). Next we introduce the span of A harmonic functions given by 

W„(f2n uj*) = span{wi £ H A (Q,r\uj*) : tOj = v t , on dui* n O U 80, n uj*, i = 1, . ..n}. (3.32) 
For future reference we note the decomposition of i? 1 (Sl n uj*) given by 

ff'(nnu*) = fl^(flnw*)+ffj(flnu*)+i 



n 



Here the first two subspaces are orthogonal with respect to the energy inner product (•, Of^nnw*)- 
The orthogonal projection from n uj*) onto H A (n Huj*) is denoted by V A . It is easily verified 

that w n {n n u*) = v A s n (n n w*). 

We now L 2 project this space onto H A (Cl n w*). The projection operator mapping L 2 (£l flw*) 
onto H A (f2 n w*) is denoted by "Po and 

ll« - ^(HMwnn) = _inf ||u - w|| L 2 (w , n0) and 

\\ v \\h(u~nn) = H^ou||i2( w *nf2) + IK 7 -^o)u||i2 (w , nn) (3.33) 

In what follows the local approximations will be chosen from the local function space VoW n (u)* n f2) 
restricted to the set w. 

As before we suppose that the side length of uj* is given by a* — a(l + p) and uj = C D £1 
where C is the concentric sub-cube of uj* and the side length of C is a. We let uj** — uj* n il 
and ct/c/2 = dist(duj n O, <9w** fl O) > 0. From the smoothness assumption on 9f2 there exists a 

dimorphism x = \I/(y) of class C 1 for which 1 S' _1 maps a/* onto uj** with <9w** n dft being a part of 
the plane y n = 0. We extend a;** across y n = by reflection. We denote the image of this extension 
under & by uj* e * and define uj* e to be given by the union uj** U (duj** n c?S7) U 

In what follows we extend u — VqV across the boundary dVL as an A-Harmonic function over uj* e 
and apply the Caciappoli inequality to recover the following theorem. 

Lemma 3.3. L 2 projection of local fields at the boundary. 

Suppose uj is a truncated cube with part of its boundary given by d£l and suppose there exists a 
dimorphism x = ^(y) of class C 1 for which maps uj** onto uj** with du>** (1 dfl being a part 
of the plane y n — 0. Let u € H A fi(Q H uj*), then there is a constant C > depending only on dfl 

such that given v € 1?(uj* n O) the projection Vqv onto H° A (f2 (1 uj*) satisfies 

II « - Po« |U (W )< ^(Z?) 1 / 2 — || u- v \\ LHu * n ci) • (3-34) 



Proof. We set u* = Pqv and and observe from Theorem 3.9 that d v (u — v*) vanishes on Oil Duj*. 
Let u = u(\&(y)) and v* — v*(\I>(y)) and note that (u — v*) is an A— harmonic function in uj** 
where A(y) = [V 1 4']~ 1 A(\I , (y))[V 1 4']~ 1 and h ■ A(u ~ v*) — 0. Here n is the unit normal to y n = 
pointing into y n < 0. Since h ■ A{u — v*) = we apply standard arguments to extend A across 
y n = so that (u — v*) is extended across y n = outside £1 as an A— harmonic function. We set 
y' = (yi, . . . ,y n -i) and write y = {y',y n )- For y n < we extend A across y n = into y n < 
according to: 1) A lJ (y',~y n ), for all i = ^', j = 1, ...,n, 2) Aij(y',-y n ) for all j ^ n,i ^ n, 3) 

—Aij(y', —y„) for j ~ n, i < n, and 4) —Aij(y' , — y„) for i = n and all j < n. 

We map back to obtain an extensio n of A across dQ and recover an A-harmonic extension of the 
function u - 



on uj e . From Theorem 



Theorem 3.1 to uj C uj* e , to find that 



3.9 



it follows that u — v* G L 2 (uj* e ) n H} oc (uj e ) and we apply 



|£(w) 



<4(/3) 



1/2 



1 



II 



V 



ap 

IU 2 (w*na) 



li 2 (a 



<4C(/3) 1/2 — || u- 

= 4C(/3) 1 / 2 — || P (w - «) \\m^nn)< 4C(/3) 1 / 2 — || u - v \\ L ^. nn) 
ap ap ' 



(3.35) 



and the theorem is proved. We point out that the analogous theorem holds for the Dirichlet problem 
and can be proved using similar arguments. 

In what follows it is always assumed that 90nw* can be flattened according to the hypothesis of 
Lemma 3.3 Next we introduce the approximation theorem associated with the space VoW n {£lP\uj*) 
given by 
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Lemma 3.4. Let u G H\ Q (fl n u>*) then there exists a v u G VoW n (Cl H to*) such that 

l/d 

||«-«u|Ua(nnu)*) = i n f II " - « IU 2 (nn«*)< ^t^"" 172 II u lU(nnw-) ( 3 - 36 ) 
v£VoW„(nnu>*) v47r 



where a* is the side length of the cube oo* , 7^ is i/ie volume of the unit ball in M. d and C n — n + 
o(l)), for d = 2, 3. 

Proof. The theorem follows immediately from an upper bound on the quotient 

it = sup mi — j— j| -. (3.37) 

Fix u G H% (Q n (J*) and for every iu G T^W^n n w*) one has g G W n (0 n (J*) such that w = Vog 
and 

„ Jk-HU^finu*) = inf J|^o(«-5)IU a (nnw-) 

ioe-p w„(nnw*) v ' gew„(nnw*) v ' 

< inf \\u- g\\ L 2/ nnu ,y (3.38) 
gew ra (nna;*) 

Thus 

D / . , I" - 5llL 2 (fin^*) ,„ „„x 

it < sup mi — jj— r -. (3.39) 

ueff^, (nno;*) sew„(nnw) ||«||£(nnu*) 

Denote the projection of u onto with respect to the energy norm || ■ Hsfnnw*) by V £ u. 

Choosing g = V e u and noting that ||(i - V e )u\\ £ (nr\u*) < \\u\\ £ ^anu*) S ives the upper bound 

d ^ IMU^nnw*) ,„ . n v 

it < sup — -. (3.40) 

ueH° A (nnw*)± w n (nna»*) IMIe(«nw*) 
Now i^ i0 (n n a;*) C #£(fi n w*) so 

D ^ IMlL 2 (nn^*) ,„ ... 

R < sup -jj — r. — -. (3-41) 

Since W„(Q ("I u*) = V A S n (fL flw*) it follows that 

{u e H^nnu*) ±r A s n (nnuj*)} = {u6^(flnu*)i5 n (finu')} 

Cjueff^n w*) _L (5„(fi n w*) + R)}. (3.42) 



Here on the second line of (3.42) the _L is with respect to the L 2 (Q, n w*) inner product. It now 
follows that 

it < Sup -jj— jj — 

, IMU 2 (nrto*) 1 ,„ .„s 

< sup irii = ; ( 3 - 43 ) 

ueH'(nnw')i (S„(ftnw*)+R) ll u ll£(onw*) v^ i "+ 1 

where /i rl +i is the largest Neumann eigenvalue associated with S n +i(u)*). One has the elementary 
lower bound fi n +i > a^n+i, where v n +i is the associated Neumann eigenvalue for the Laplacian on 
bj* n f2. For this case Weyl's theorem [28] gives 

= 4tt f " ) + o(n" d ). (3.44) 
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The upper bound on R now follows from (3.441 together with the inequality |fi fl u>*\ < (a*) d and 
the theorem is proved. 

We introduce the local space near the boundary given by dF^uj* Dfl) = VoWr^ui* n O) and 
define the local approximation space dJ 7 n (uj,uj* fl O) to be given by the restriction of d!F m {uj* n fi) 
on uj. 



Now we apply Lemma 3.3 to u — v u on u) C l~l u* and combine it with Lemma 3.4 to obtain 
the following convergence rate associated with the family of approximation spaces dJ- n (uj, uj* (1 57). 

Lemma 3.5. Let u G H^ (uj* fl £7), then there exists an approximation v u G dJ- n (uj,uj* n f2) /or 
which 

\\u-v u \\s(u)= ini II tt-w |U( W ) < nfi)C n || u Hs^.nn) (3-45) 

w/iere 

I(w,u>* HO) = 8C^V ^^(/3/a) 1/2 and C„ = n" 1 ^^ + o(l)), d = 2, 3, (3.46) 

where is the volume of the unit ball in dimension d and C depends only upon dQ. 

Now we proceed iteratively to construct a family of local approximation spaces with a rate of 
convergence that is nearly exponential. For any pair of two concentric cubes Q C Q such that their 
intersections ui = Q n 57 and uj = Q n 57 have nonzero volume we define dT n (u), cD) to be the space 
given by the restriction of VoW n (uj) on u>. We recall that the two concentric cubes C C uj* are 
of side length a and cr* = er(l + p) respectively and ui = C fl Q, see Figure [3] Let iV > 1 be an 
integer and consider the nested family of concentric cubes Qj, j = 1,2, ...N + 1 with Qj+i C Qj 
and Qi = and Qjv + i = C. The side lengths of Qj are given by ct(1 + p(N + 1 - j)/N))/2. Set 

= Qj H 57 to obtain w = w^v+i C lo^ C ■ ■ • c ui = n 57. We introduce the local spaces, 
dJ r n (uj,uiN), dJ 7 n (uj,uiN-i),---,dJ 7 n (u),u!i). Put m = N x n and we define the approximation space 

* = ar(m,w,w*nl]) =aF n (w,u; 1 ) + -"+aF n (>j,u> A r)- (3.47) 

The convergence rate associated with the local approximation space dT(m,oj,uj* Hfi) is given in the 
following theorem. 

Theorem 3.8. Let u G Hl (u>* n £7), then th ere exists z„£$ = dT(m, uj, uj* (1 57) smc/i i/ia£ 

||«-*«. || f(w) <^ || u |U (w .nn) (3-48) 

and ? = 8C ^±±£jV(/3/a) 1/2 C- n . 

Proof. The proof is by induction and is identical to the proof of Theorem |3.5| 
Theorem |3 . 7| now follows on choosing the appropriate N and using arguments identical to those 
used to establish Theorem 13.31 

We conclude by stating and proving the following theorem. 

Theorem 3.9. The set H A (u>* flft) is a subspace of the space of A-harmonic functions belonging 

to L 2 (fl flw'). Functions v belonging to H a,q(ui* H 57) have the following local properties. For any 
open subset O d uj* H Q, 

1. V £ H l {0) for any O C uj* n £7, such that dist(dO, du* nU)>0, and 

2. ifdOn(dnnuj*) ^ then d v v = on dO n (<957 n w*). 

Proof. Given G o( w * ^ ^0 then there is a sequence it n G fl^ o( w * H O) such that u„ — > 
in I 2 (w*flf!). We show first that Uoo is ^4-harmonic on u'flO and belongs to H} oc (uf fl£7) . To see this 
pick any ball B{x$,r) CC a/ n 57 centered at of radius r. We apply the Cacciappoli inequality 
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(Theorem 3.1 1 together with the Rellich-Kondrachov compactness theorem to deduce that u n is 
Cauchy with respect to the energy norm in B(xo, r/2). From the completeness of H 1 (B(xq, r/2) we 
see that u n — > Uoo in H 1 (B(xo,r/2)). From this we conclude that Moo € H} oc {w* n 51). The weak 
formulation of the boundary value problem together with the strong convergence of the sequence 
easily shows that is A-harmonic. 

Next consider any open subset O C cj*n£l such that dist(dO, <9w*n51) > and dOn(dQnu)*) ^ 0. 
Consider any ball B(xq,t) centered at xq £ 951 of radius r with B(xo,r) n 51 contained inside O. 

The dimorphism x = ^~ 1 (y) maps B(xo,r) fl 51 onto B(xo,r) H 51, with B(xo,r) H 951 being part 

of the plane y n = 0. Extend i?(xo,r) n 51 across y n = by reflection. Denote the image of this 
extension under * by C* E and define C* = (B(x Q , r) n 51) U C^. Now consider « n £ if^ (w* n 51) 

for which u n — > in L 2 (u>* PI 51). Now u„ can be extended as an A— harmonic function over C* 
with || tin. 1 1 i 2 (c* ) < C||' u ri||L 2 (s(2;o.r)na) where C depends only on 951. Since B(xQ,r/2) n 51 C C* 



we can apply a Cacciappoli inequality analogous to Theorem 3.1 to discover that {un}$?Li is a 
Cauchy sequence in H 1 (B(xq, r/2) n 51) and we conclude that u n — > Uoo in H 1 (B(xo, r/2) nil). This 
establishes property (1) of the theorem. 

Observe that since 9„it n = n ■ AS/u n vanishes on B(xq, r) n 951 we can write 

\\9uU o\\H-' L / 2 (B(xo,r/2)n8n) = II ~ C?i/«oo || H- 1 / 2 (B(x ,r/2)ndfl) 

< C\\u n - W 00 ||fl-i(B(xo,r/2)nf2) (3.49) 

where C is independent of n. Property (2) now follows on noting that 

lim \\u n - Uoo||ffi(B(x ,r/2)nn) = 0. (3.50) 



3.3 Nearly exponential convergence for GFEM applied to heterogeneous 
systems 

Theorems |3.3| and |3.7| provide the local finite dimensional subspaces required for a global Galerkin 
approximation with error that converges nearly exponentially with the degrees of freedom. For a 
given partition Wj,t = l,...,mwe denote these subspaces by i = 1, . . . , m. We augment each of 
the subspaces with the subspace of constant functions denoted by K and write VE^ = $j © R. For 
domains that touch the boundary of 51 the local approximations are taken from the hyperplane 
© u^. Here u p £ Ha(oj* D 51) is the local particular solution introduced in the previous section 
that satisfies n-AV u v = g on w*n951. We denote the dimensions of by JVj and set N = Y^Lx N%- 
Recalling Theorem |2.1| and applying Theorems |3.3| and |3.7| we obtain the following approximation 
theorem. 

Theorem 3.10. Nearly exponential approximation for GFEM. 

For e > there is a N E > such that for all N > N E , there exist Q £ ^ for uji n 951 = 0, 
Ci G © Wp, /or Wi H 951 7^ and a constant K, independent of N such that the approximation 
Q £ _ff 1 (51) given by 



satisfies 



and 



II Uq - C IU 2 (n) < £' 



_Ar(i+d e ) 



(3.52) 



l|w -CIU(n) < fa 



_Pf{i+d e ) 



(3.53) 
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4 Implementation of the multiscale GFEM method 



In this section we provide an overview of the main ideas noting that the specific challenges and 
details of the implementation are the focus of future work. The implementation consists of three 
parts: 

• Multiple independent parallel computations for construction of the local bases and the subse- 
quent assembly of the global stiffness matrix. 

• A single global computation using the global stiffness matrix and load vector. 

• Recovery of preselected local features of the solution through the multiplication of the local 
bases by solutions of the global problem, e.g., the recovery of stresses at fiber matrix interfaces. 

We now give an outline of the primary issues involved in the computation of the local optimal 



approximation spaces ^F;, i = 1, • • • , m prov ided by Th eorems 3.3 and 3.7 The local bases are given 



by the eigenf unctions of the problems (3.6 1 and (3.301. In what follows we will assume that all the 
sub-domains uji are roughly the same size and we will suppress the index i and write lo and \£ . 

A suitable and effective numerical method for the construction of the local basis '5 is given by the 
subspace approach (see [21] chapter 11). This method is based on the Raleigh-Ritz approximation. 
Here the key ingredient for the success of this method is the selection of a suitable subspace with 
span that should approximate the span of the first iVj eigenfunctions. We now briefly discuss the 
construction of the subspace and the discrete representation of the eigenfunctions used in the local 
basis of dimension N > 1. We start by introducing M > N functions ft , k = 1,2, .. defined 
on the boundary dto* . The example given at the end of this section shows that good candidates 
for ft are the normal derivatives of the harmonic polynomials of degree k. Other choices are also 
proposed in section 6. These functions are then used to construct the M dimensional subspace & of 
A-harmonic functions u^ k on w* which satisfy the Neumann boundary condition n ■ AS/u = ft. This 
subspace © is used within the subspace approach to construct the desired N eigenfunctions and 
eigenvalues. These eigenfunctions will comprise the local approximation ^ used in the multiscale 
GFEM. The appropriate selection of N and M is determined by the rate of decay of the eigenvalues 
with respect to these parameters. The numerical construction of the local basis for each subdomain 
can be carried out in parallel. 

Since we have established nearly exponential convergence it is expected that only a small number 
of eigenfunctions will need to be computed. Moreover the examples presented in subsequent sections 
show that the desired eigenfunctions can be smooth on duj* allowing for an accurate approximation 
for relatively small values of M > N. 

The global basis is constructed by combining the local bases with the partition of unity functions. 
The partition of unity structure guarantees a sparse global stiffness matrix and the assembly of the 
global stiffness matrix is also local procedure that can be carried out in parallel. 

We now provide a rough estimate for the computational work involved in the GFEM for het- 
erogeneous systems for problems posed over a computational domain Q C R 2 . To start we cover 
to* by a finite element mesh with elements of size h and solve the Neumann problem for the differ- 
ential equation div(AVu) — 0, subject to the boundary condition ft. These solutions deliver the 
(approximate) A-harmonic functions vfi. For example using Gaussian elimination (LU decomposi- 
tion) we require (| uj* \ hr 2 )~ 2 operations (due to the sparseness of the stiffness matrix). Since 
M << h~ 2 the cost of the computing M functions does not change the order of operations. Fur- 
thermore computation of the N eigenfunctions is relatively small and hence the cost of the creation 
of the space ^ is of order | |~ 2 h~ 4 . The computation of the entries in the associated stiffness 
matrix will not change the order of operations. The major problem is the choice of mesh size h that 
leads to an acceptable accuracy. If the boundary functions ft are smooth, for example the traces 
of harmonic polynomials as mentioned above, then we conjecture that the accuracy is on the order 
of e a =|| w Sfc — u^ k \\s(u)< Ch 1 M s || u In future work we plan to analyze this conjecture. 

For the case of a proper mesh applied to the fiber material discussed in section 2, we expect 7 = 2 
and the effect of M is negligible. The span of the functions it Sfc are then used to construct the 
approximate space ^> using the subspace method. These approximate spaces are then used within 
the GFEM scheme. 

On applying Gaussian elimination to the global stiffness matrix we obtain the approximate 
solution over fl. The local error incurred by approximating the solution u over a subdomain uj using 
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the approximate local basis ^ h will be of order e™^ 1+d . The size of the global stiffness matrix 
is Nx | D, | / | u \. Because the stiffness matrix is sparse, the work of solving the global system 
using Gaussian elimination is (Nx | f2 | / | u> |)~ 4 . The final implementation issue involves the "best 
choice" of the size of the local domains lo and lo* for maximum computational efficiency within the 
context of parallel computation. This question is not addressed here however it is clear that it is a 
very important issue that needs to be addressed in the implementation. 

We conclude with an example that illustrates the exponential convergence. In (26] the numerical 
solution for the deformation inside a shaft reinforced with long compliant fibers with zero rigidity is 
given. The material in between the fibers is referred to as the matrix. Here the shaft is subjected 
to anti-pane shear loading and the system of elasticity reduces to the single scalar equation for 
deformations u perpendicular to the mid-plane of the shaft. When no fibers are in contact with each 
other the entire theory developed here also implies the exponential convergence of the GFEM for 
this problem. The computational domain is the shaft mid-plane f2 given by a subset of the x — y 
plane portrayed in Figure |4| For this problem when one constructs the local basis over a generic lo* 
the convention is to remove any fiber domains intersecting duo* and to replace with matrix material. 
For this kind of problem the local basis functions u ffc are taken to be harmonic in the matrix outside 
the fibers, taking zero Neumann data on the boundary of the fibers and taking Neumann data on the 
boundary of duj* given by the traces of harmonic polynomials of degree k. The Neumann condition 
du/dn = g posed on the boundary of the computational domain dfl is given by g = (2x — y). 
We now describe the lo.i comprising the partition of unity for this example. The computational 
domain Q is covered by a 16 x 16 mesh of square elements r. The partition of unity functions are the 
standard "hat" functions of the finite elements which are bilinear on the elements r. The supports 
of these "hat" functions create the set of local domains uii. Sets LOi interior to f2 are composed of 
2x2 squares. The interior domains lo* are composed of 6 x 6 squares. The domains LOi and lo* 
close to the boundary of are constructed as described above in section 3.2. The spaces Vf^ are the 
restrictions to loi C lo* of the solutions M ft on lo* with Neumann boundary condition given by the 
traces of the harmonic polynomials of degree k, hence the dimension of 'Pj is 2k + 1. These local 
basis functions are computed using the finite clement method on sufficiently fine mesh, so that the 
error is negligible. 

The "exact" solution u is computed by "overkill." The relative energy norm of the error as 
function of the degree k of the harmonic polynomials is presented in Table 1. Figure [5] shows the 
error in the log scale as function of k in the linear scale. The straight line clearly shows an exponential 
rate of decay. From these numbers we see that the rate is given by e -° ASn while the estimate is 
given by e -0 - 33 ", where n = 2k + 1. 
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Figure 4: Mid-plane of shaft. Reprinted with permission from [2B]. Copyright 2004, John Wiley and 
Sons. 

The simulations presented in [26 also numerically investigate the effect of the distance between 
the boundaries of lo and lo* . There it is found that the rate of convergence decreases as the distance 
is reduced and that the exponential convergence vanishes when lo and lo* coincide. 
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0.89% 
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Figure 5: Decay of the error. 



Table 1. The error as function of the degree of harmonic polynomials characterizing the boundary 

conditions q k . 



So far we addressed only the case that the right hand side / = in the equation (1.1). The 
general case / 7^ can be easily reduced t o th e case / = 0. To this end we introduce the local partic- 



ular solution of the differential equation (1.1 ) on each ui* subject to a constant Neumann boundary 
condition over dw* determined according to the consistency condition. This is inexpensive to imple- 
ment because the stiffness matrix has already been constructed together with its LU decomposition. 
We denote the local particular solution by u Ui . Then the local approximation over uii used in the 

GFEM belongs to the hyperplane © u uli . Here the finite dimensional subspace is given by the 
optimal local bases constructed for the A— harmonic problem. As before this construction delivers 
a nearly exponential rate of convergence. 

5 Homogenization of the n-width and exponential decay of 
approximation error in the pre-asymptotic regime 

We identify the homogenization limit of the n-width and the corresponding optimal basis functions. 
These ideas are used to provide examples of exponential convergence of the approximation error 
when the characteristic length scale describing a heterogeneous medium is sufficiently small. In 
what follows we work in the general context and homogenization is described by if-convergence 
[17] or G-convergence [23] ■ We consider a sequence of coefficient matrices A e (x) in £ indexed by 
e, with e = l/£ for £ — 1,2,.... Since we consider symmetric coefficient matrices the notions of 
G convergence and H convergence coincide and the class of coefficients £ is compact with respect 
to iJ-convergence see P~7J, [33] ■ In what follows we assume that the sequence A e if -converges to a 

homogenized coefficient matrix A° in € and we write A e — > A . 

We describe the n- widths associated with the sequence A e and the if- limit A . For each A e we 
introduce the Hilbert space H t defined to be all elements in Ha'{uj*)/M. equipped with the energy 
inner product 

(u,v) Hc = J A c Vu-Vvdx (5.1) 

and norm \\v\\^ = (f,f)« e . The Hilbert space associated with ^4° is denoted by Ho and is defined 
to be all elements in H^o (w*)/K equipped with the energy inner product 

(u,v)n = A°Vu-Vvdx (5.2) 



and norm ||w||^ o = {v,v) Ho . 

For each A e we introduce the restriction operator P e : H e — > Ha* (to) /K such that P € (u)(x) = u(x) 
for all x £ uj and u £ T-i t . As mentioned in the previous section the operator P e associated with A e 
is compact. For future reference the energy bilinear form defined on H^e(a;)/K is given by 

(u,v) Ec ^= / A e Wu-Wvdx (5.3) 

J UJ 
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and we set ||tt||| = (u, w)f e (u)- Similarly, for the P-limit A we introduce the compact operator 
given by the restriction P : Ho — > Hao(co)/M. such that Pq(u)(x) = u(x) for all x € uj and u G Ho- 
The energy bilinear form defined on H^o (ui) is given by 

(u,v) £o{u) = [ A a \7u-\7vdx (5.4) 

and we set ||u||| o(w) = (u,u) £o(u) . 

The n width associated with the coefficients A e is given by 

£(w,w*) = inf sup inf H P ^~*I^(") . (5 . 5) 

S(n)CH At (u,)/K uen€ xeS(n) \\u\\ Ht 

The optimal local approximation space ^(cj) associated with A e is described in terms of the 
cigenfunctions associated with the following spectral problem. We introduce the adjoint opera- 
tor P* : Ha'{u>)/M. — > H e and the operator P*P e is a self adjoint non-negative compact map taking 
H e into itself. The eigenfunctions and eigenvalues are denoted by {<-p\} c £L 1 and {Af}-^ and satisfy 
the problem 

(P;P e ^,<5) Ke = A?(tf,<J) We ,We« e . (5.6) 

The non-zero eigenvalues of P* P e are listed according to decreasing order of magnitude 

X{ > Xi, > . . . > 0. (5.7) 

The optimal approximation space is given by ^(w) = span{ipl, . . . ,ip„}, where — P e (fj and 
d^(ui,uj*) = ^X f n+1 . The n width associated with the coefficient A is given by 

jO/ *\ • r ■ r \\Po U ~ Xl Uo(w) /- „n 

<ft[uj,(jj )— mf SUP inf rr— T . — . 5.8 

S(n)CH A o(u) ue -H oX eS(n) \\u\\ Ho 

The optimal local approximation space ^(w) associated with A is described in terms of the 
cigenfunctions associated with the following spectral problem. We introduce the adjoint opera- 
tor Pq : Ha' (w)/M — > Hq and the operator Pq P is a self adjoint non-negative compact map taking 
Ho into itself. The eigenfunctions and eigenvalues are denoted by {tp®}^ and {X®}°1 1 and satisfy 
the problem 

(P*P ^,8) Ho = X^lS^ySeHo. (5.9) 

The non-zero eigenvalues of P *Po are listed according to decreasing order of magnitude 

A?>A°>...>0. (5.10) 

The optimal approximation space is given by ^(w) = span^i, . . . , ^°}, where ip® = Pofi and 
d°(cj,cj*) = y/ A° +1 . The homogenization limit of n- widths and optimal approximations is given by 
the following theorem. 

Theorem 5.1. Suppose that the coefficient matrices A e {x) in € H-converge to A°(x) in <£ as e — > 0. 
TTien i/iere exists a subsequence of coefficients A e -5- A suc/i t/iat 

A?-S-A° and ip\ -± ip° t , weakly in H 1 ^*), for i = 1,2 .. . (5.11) 



-Hence 

lim(£(a;,w*) = (5.12) 



e->0 



and eac/i function in the optimal basis for A € given by ^ (to) = span{tp\ , . . . , tf)^} converges weakly in 
H 1 ^) to the corresponding function in the optimal basis for A given by (w) = span{ip®, . . . , -0°}- 
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Proof. The proof proceeds in three steps. 

Step 1. We start by fixing the index i and state the following Lemma. 

Lemma 5.1. Suppose A e A° . For i fixed, consider the associated eigenfunction and eigenvalue 
(ipj,Xl) of P*P e . Sending e — > and passing to subsequences as necessary there exists a positive 
number Xi and function Tp i £ Hq for which 

XI — > Xi and ip\ — ^ weakly in H 1 (cj*), (5.13) 

and 

(P*P ^,6) no = \i(Vi,6)H , VSeHo- (5.14) 



Before proceeding with the proof of Lemma |5.1| we state the following two compensated com- 
pactness results presented in the work of Murat and Tartar [T7] for later reference. 

Lemma 5.2. Let e — >■ 0, D be any open subset of Mr, and £ £ , v e , be sequences such that 



C G L 2 (D) d , 

f £° weakly in L 2 (D) d , 

div^ e — > div^° strongly in H^ 1 (D), 



(5.15) 



Then 



v e — 1 v° weakly in H 1 (D). 



[ C-^v^dx^ [ £° • \7v°ndx, V77 € C$°{D). 
Jd Jd 



(5.16) 



(5.17) 



Lemma 5.3. Suppose that A e (x) in £ H-converges to A°(x) in £ as e — > 0. Assume that u e € 
if *(£)), f € H- X (D) and 



dw(A e Vu e ) = / e , in A 

u £ — u° weakly in H 1 (D), 
f -> /° strongly in H~ 1 (D), 



(5.18) 



/or e — > 0. TTierc 



/ (A e Vu e • Vu^ndx 
Jd 



A e Vu e ■ Vu e 



(A°Vw° ■ Vu Q )ndx, Vry <= Cg°(D) and 
^°Vm° • Vw° weafcfy in L} oc {D). 



Proof of Lemma 5.1 Following section 3 the we write (5.6) as 



A e Vipt -V5dx = Xt I A e Vip\ ■ V<5 dx, V<5 e % 6 



(5.19) 



(5.20) 



Now consider the sequence of eigenfunctions {ifl} e> o for (5.20) and with out loss of generality 
we normalize <p\ so that 



A^ipl ■ \7ip\ dx = 1 and Af = / A e V^| • V<^ dx 



(5.21) 
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From (5.21| we can extract a subsequence {v!}e>o and (p i £ i? 1 (o;*) such that ipj ip i weakly in 
Since tp% £ % e C we apply Lemma 5.3 to deduce that for e — > 0, 

A £ V<^ £ A°V^ weakly in L 2 (uj*) d , 

I {A e \7ipl-\I^)8dx ->> I (APVipi'VvJSdxW eC °°(w*) and 

A £ V^i-Vv- A°V^i • Vp <; weakly in L 1 ^), (5.22) 

and <^ G Hq, hence 

Af -> / A°V<^ • V<^ = A*. (5.23) 

To finish the proof we show that ^ and Ai are solutions of ( 5. 14 1 . Consider any g £ H 1 ^ 2 (dui*) 
and (S c € "H e such that <5 £ = g on We write S e = w e + v where 

v = g on did* and Aw = in lu* , (5-24) 

and w £ € i?o(u;*), where 

- div(A £ Vw £ ) = div(A £ V«), (5.25) 



For any sequence of coefficients A e E £ Theorem 1 of [TB] shows that the sequence Vu> £ enjoys 
the higher integrability given by the following Lemma. 



Lemma 5.4. There is an interval Q > p > 2 such that for v £ W 1,p {lu*) ther 

sup{||Viy e || £ , P ( a ,.)} < oo. 

.ffere £/ie interval is independent of e and depends only on lu* , a and /3. 



(5.26) 



Now consider a dense subset S of H x ^ 2 (dLU*) such that g £ S implies that the sol utio n v of (5.24) 
belongs to W l ' p (w*). Then ||V(5 e ||ij.( w ») < ||Vw £ || L p( u .) + ||Vu||lp( w ») and Lemma 5.4 implies that 
the associated sequence {<5 £ } c >o satisfies 



sup{||V<5 £ || L p (u ,)} < oo. 

e>0 



(5.27) 



Additionally passing to a subsequence if necessary we see that the re is an element S £ _ff 1 (w*) for 
which S e — 1 6 we akly in H (ui*). Next an application of Lemma 5.3 shows that 5 £ T-Lq and an 
application of 5.2 to the sequences {A £ V Vf} e >o and {S e } e> o gives 



/ {A^Vipl ■ \7S e )ndx -> I (A°V<p i -VS)r)dx'ir) £C^{u 

J LJ* J U)* 



From (5.27) wc deduce that {A e V(p\ • V(5 € } e> o is cquiintcgrablc on lj* and it follows that 



/ 



A'Vvl ■ V£ £ 
A e Vipt • V<5 £ 



From (|5.23|), (|5.29|, (|5.30|), we deduce that 

A°VTp i -V5dx = A, 



/ A°VXpi ■ \7Sdx, 

J LU 

[ A°S7^ r \75dx. 

J LU* 

A^lfii-WSdx, 



(5.28) 

(5.29) 
(5.30) 

(5.31) 
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or equivalcntly 

(P *Po^,<5) Mo - kfaiQno (5-32) 



for all test functions S belonging to "Ho with traces in S, i.e., S = g on duo* , for g 6 S. Lemma |5.1 
now follows from the density of 5 in H 1 ' 2 (doo*) 



Step 2. We apply Lemma 5.1 together with a diagonalization argument to extract a subsequence 



still denoted by (Af , tpf), such that for every i = 1,2. 

A- -> Ai and ^ ->> ^ weakly in H^uo*) (5.33) 



where (Xi,(pA are solutions of (5.14) and 

Ai > A 2 > . . . > 0. (5.34) 

Step 3. The final step is to show that all the eigenfunctions and eigenvalues of the operator Pq P 
are given by (A,,^) obtained in step 2. We argue by contradiction and assume that there is an 

eigenvalue A of Pq Pq for which A 7^ Ai for every i = 1,2, Let up be a corresponding normalized 

eigenvector i.e., (Pq PqP, §)n = A(<,o, <5)-h , for every S £ Hq and ||y||w = 1- Then there is an 
integer m such that 

A > A m+ i. (5.35) 
To proceed we introduce the Rayleigh quotient for v g H e given by 

(P; P e v,v) n. 

and the eigenvalues of P*P e listed in decreasing order are given by 



Re{v) = ^i^^ (5.36) 



Af = max R e (v). (5.37) 

»EW t lpi,...?!|_ 1 

We establish the contradiction first under the extra assumption that the gradient of <p € Ho enjoys 
higher integrability and belongs to W 1,p (uo*) for p > 2. We then indicate how to proceed without 
this assumption. 

Introduce u e e "H e such that u e = up on On passing to a further subsequence if needed we 

apply Theorem |5.3| to see that 

u e — ip weakly in H l (io*) and 

(P*P e u e ,u € )n t = / A £ Vi/ • Vu e ->■ A A°\7tp ■ \7ipdx = A. (5.38) 



Noting from that y> G VK 1,p (w*) we observe from the arguments preceding (5.27) that 

sup{||Vu e || LP(w ,)} < 00. (5.39) 

Thus the sequence {A £ Vit £ • Vu e } e> o is equiintegrable on a;* and we conclude that 

(u £ ,M £ )-H e = / A £ Vu £ ■\7u e dx -> [ A°S7ip- \Jtpdx = 1 (5.40) 

J UJ* J UJ* 



lfmP e (u £ ) = A. (5.41) 

Now introduce t> £ € "H e given by 

m 

»=i 
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As before we make use of the equiintegrability of {A e Vu c ■ Vp\} t> o on w* together with Lemma 5.3 
to find that 



I A £ Wu e ■ Vyj| dx-> A°Wip ■ Vifi dx = (ip, p,)^, 

J UJ* J CO* 



(5.43) 



Since A 7^ A, for all i, it follows from the orthogonality of eigenvectors of PqPo that (ip,ip i )u = 
for i = 1, 2, . . . , m and we deduce that 



(u e ,<pl) n€ ->0. 



Writing 



(5.44) 



ll«1«. = ll«%.-£( Ue >tf 

m 

(P:P e v\v^) nt = {P*P e u e ,u e )- H<i -^\l{u € ,ipf)\ 



i=l 



and sending e to zero using (5.38), (5.40), and (5.44) we conclude that 



On the other hand 



lim R e (v e ) = A. 



(u £ ,<^)« e =0, fori = 1,2, 



(5.45) 



(5.46) 



(5.47) 



so from (5.37) we get A^ +1 > A and taking limits gives A m+ i > A which is a contradiction to the 
original assumption A > A m+ i. 

We now remove the higher integrability assumption on the gradient of <p € "Ho- For this case 
consider a sequence s — 1/1,1= 1,2,... and functions 6 S € W 1,p (oj*) that converge to ip in W 1,2 (uj*) 
as s goes to zero. Choose u\ € H e such that = <5 S on duj* . Then construct according to 



v: = u: 



(5.48) 



As before (wj, (pf)n t = 0, for i = 1, . . . , m and A^ n+1 > i£ e (uj). Following previous arguments one 
deduces that the sequence u £ s is bounded in W 1,p (lo*) and on passing to subsequences as necessary 

u e s u s weakly in H (oj*) where u s G TL , 

(P* t P e ul,ul) Hc ^ f A°Vu s -Vu s dx, 

(«>Dh« -> (u s) U s ) no , and 

{u e s ,vt)-Ht -> (u s ,<Pi)u , for i = l,...,m 



and 



lim i? e (vj) 



(5.49) 



(5.50) 



Since 5 S converges strongly in iJ 1 (a;*) to ^ it follows from the uniqueness of solution of the Dirichlet 
boundary value problem for A harmonic functions that u s converges strongly to tp in TLq thus 



A m +i > lim lim R e (v e ) = A 

s->0e->0 



(5.51) 



and we arrive at a contradiction and Theorem |5.1| is proved 
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We conclude by applying the homogenization of n-width theorem to construct an example that 
shows exponential decay of the approximation error in the pre-asymptotic regime. We consider a 
heterogeneous medium with characteristic length scale e > 0. To fix ideas we work in two dimensions 
and suppose that the associated sequence of coefficients A e is such that it fZ-converges to a constant 
effective conductivity A° matrix as e — > 0. In the coordinate system corresponding to the eigenvectors 
e 1 , e 2 of A we have A = axe 1 <£> e 1 + a^e 2 ® e 2 and we set b = 0,2/0,1. To fix ideas we suppose 
that uj* is the ellipsoid E r * = {(xi, £2); x\ + x\/b = r*} and uj C uj* is the concentric ellipsoid 
E r = {(ii, X2); x\ + x\/b = r} with r < r*. For z = x + iy recall the harmonic polynomials 
Wj(xi, X2) — ytzi — ri cos (j0), Wj(xi, x%) — ^$z 3 — r J sin (jO), for j = 1, . . . , n. Calculation shows 
that the optimal basis associated with the n width for A is given by the ^4° harmonic polynomials 
vj = Wj(xi,X2/Vb), Vj = Wj(xi, X2/ Vb) and eigenvalues Xj — e -2 ' 1 "^^, j = l,...,n of 

/ A°V(pj -\75dx = Xj [ A°Vcpj -\75dx, (5.52) 



for all 5 G 7/40(0;*). It follows from Theorem 3.1 that the decay of approximation error for the 
optimal basis associated with the homogenized coefficient A is 



(5.53) 



Now we denote t he n width associat ed wi th the optimal basis for A e by d e n (E r , E r »). Direct applica- 



tion of Theorem 5.1 together with (5.531 gives the the following bound on the pre-asymptotic rate 



of approximation error. 

Theorem 5.2. Given N > and tolerance r > there exist an e > such that for 1 < n < N, that 

e -|ln^|(„+l) _ T < d ^ Er}Er ,) < e -\hx^\(n+l) +T _ ( 5 54 ) 



6 Implementation in the pre-asymptotic regime and more 
examples of exponential convergence 

In this section we discuss a method for computational approximation that employs the optimal basis 
for the homogenized problem to construct approximation spaces for composites with heterogeneities 
on the length scale e > relative to the size of u>* . We work in the general context and consider a 
sequence of coefficient matrices {A e } e> o £ £ that Jf-converg e to a homogeni zed c oefficient matrix 



A £ €. For this case we recall the eigenfunctions ip\ of (5.6 1 and $ of (5.9) associated with 
A e and A respectively. For e > fixed the optimal approximation space is given by the span 
of the restriction of the functions ip\, i — 1, ... ,n to uj. However in general it is known that the 
direct numerical computation of eigenfunctions is computationally expensive. Instead we introduce 
the functions 4>\ £ Ha^(^*)/R such that </>f = (p° on duj* , for i — 1, . . . ,n. We then define the 
approximation space V e n (w) by 

V?{oj) = span{ul = Pct>t, i=l,...,n} (6.1) 

and state the following approximation theorem 

Theorem 6.1. Given a tolerance r > there exists an e > such that Ve < e 

IK - <Pi\\e(u) < T (6.2) 

We point out that this theorem remains the same if we choose (\>\ £ Ha'(oj*)/R such that 
n ■ A £ V(f)\ = n ■ A^Vtp® on dw* , for i = 1, . . . , n. When the homogenized coefficient A is sufficiently 
simple e.g., A is a constant, and oj and uj* are concentric ellipsoids, the optimal approximation 
space for the homogenized problem is given by explicit transcendental functions. And it follows 
that the associated approximation space V n (uj) is far less expensive to compute than the eigenvalue 
problem associated with the optimal approximation space. For these situations Theorem |6 . 1 1 shows 
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that V™(ijj) can be used provided that e is sufficiently small. We point out that the traces of the 
approximations </>f are indeed smooth on 8lu* noting that this is exactly the assumption made in 
section when considering the accuracy of the approximate local basis given in section [4j For fiber 
reinforced composite materials it is clear that the size of lu needs to be chosen sufficiently large so 
that the relative length scale of the fi ber c ross sections as chara cteri zed by e is sufficiently small. 

We now give the proof of Theorem 
ip\ -> <p\ in L 2 (uj*). On the other hanc 



G.i 



Recall from Theorem 5.f that p\ tp® in H 1 (lu* ) 1 h ence 
since A e _ff-converges to A" it follows from Theorem 5.3 that 
tpl in iJ 1 (o;*), hence <f>\ — ¥ ip® in L 2 (lu*). Application of the Caccioppoli inequality delivers 

IK - tflkuo < m 1/2 /* P )Hl - <pt\\» (u .) 

< mf l2 lap) (||0? - d\\L* {u *) + M- P?IU»( W .)) (6-3) 



and Theorem |6.1| is proved. 

In the numerical example presented at the end of section|4 we have assumed that the homogenized 
equation is given by the Laplace equation, i.e., A° — I and that the functions cp® = q are the traces 
of harmonic polynomials. 

Consider a family of heterogeneous media with characteristic length scale e > 0. We suppose as 
before A e is ^/-convergent and converges to a constant effective conductivity A matrix as e — » 0. 
We take lu* to be the unit square and w to be a concentric square of side length a < I contained 
inside lu* . We suppose that a is such that we can fit concentric ellipsoids E r C E r », with r < r* 
inside lu* such that u) is contained inside the smaller ellipsoid E r . We consider even dimensional 
approximation spaces and take our approximation space V e ™(w*) to be given by the span of the 
A e -harmonic functions (j)j on lu* taking the Neumann data given by n ■ A°V«j, for j = l,...n/2 

and 4>j taking the Neumann data given by n ■ A°'Vvj, for j = 1, . . . , n/2. Here Vj and Vj are the A 
harmonic polynomials introduced in the previous section and n is the outward directed unit normal 
on the boundary of lu* . For this case we have the following theorem. 

Theorem 6.2. For any sequence {w e } e >o C Ha'{lu*)/M. such that sup e> g{||u e ||£ e ( w *)} < oo, then 
given n > and tolerance r > and on passing to a subsequence if necessary there exist an eo > 
such that for e < 6q 

inf J|x-u«lkM<(e- |Jn ^ l(n+1) +r)K|U. (<l ,. ) . (6.4) 



Proof. Let {ipf, . . . , ^/ 2 , i>\, ■ ■ ■ , ^ e ^ ne optimal basis for the concentric ellipsoids E r C E r * 

for the coefficient A e . The subspace spanned by these functions is denoted by W"(E r *). The optimal 
basis for the concentric ellipsoids E r C E r * for the homogenized coefficient A , denoted by W a (E r *), 
is precisely the span of the A harmonic polynomials Vj = Wj(xi, X2/ Vb), Vj = u>j(x±,X2/Vb), 
j = 0, . . . , n/2. Then there is a sequence of constant vectors {cf , . . . , c? , 2 , cf , . . . , , 2 } bounded in 

MP such that V, G WP{E r *) is given by = E^Lli^j + and for Xe = E"=i(c$<^ + 

we deduce that 

inf ||x - u e \\ £c{u) < \\ X t - u e \\ £t{Er) 

< \\U £ - Ipeh^Er) + llXe - IpeU^Er) 

< <%{E r ,E r »)\\u e \\e € (E r .) + llXe - i>e\\e t (E r ) 

< d n e {E r ,E r ,)\\u e \\ £c{Ert) + llXe - i>e\\LHE r ,). (6.5) 

Here the last inequality in ( |6.5[ ) follows from Theorem 3.1 and S — dist(dE r *,dE r ). Moreover since 
{(cf , . . . , c^, cf , . . . , c^/ 2 } is bounded in R n we can extract a convergent subsequence and from our 

previous observations on H convergence we have that there is a ipo G WP(E r .) such that Xe - ► 4>o 
and tjje — > ipo strongly in L 2 (E r *). It now follows that 

||Xe - Uellfi.^) < dP{E r ,E r ,)\\u e \\ £rXEr , ) 

H c — (\\A ~ L 2 (E r * ) + llXe - AWl^e,,,)) , (6-6) 
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and the the orem is proved. 

Theorem 6.2 shows that the use of V™(uj*) delivers exponential convergence in the pre-asymptotic 
regime when the size and separation of the disks is sufficiently small. 




Figure 6: Level lines of fiber volume fraction. 



These examples demonstrate how homogenized coefficients can be used in the construction of 
the optimal shape functions. Of course the question of the "best" choice of homogenized coefficients 
which lead to a reasonable approximation for general situations is not clear. Nevertheless one can 
formally proceed by selecting appropriately sized uj* such that it is large with respect to the features 
of the heterogeneity but such that the heterogeneity is statistically uniform within it. With this in 
mind we return to the fiber composite portrayed in Figure [T] Figure [6] is a map of the level lines 
of the volume fraction taken over a moving window given by the square of side length 116 /zm, [2]. 
It provides a characterization of the spatial variability of the material. The volume fraction varies 
between 45% and 65% across the sample and clearly demonstrates the statistical inhomogeneity of 
the material. The correlation between volume fraction and effective elastic properties for this sample 
is illustrated in Figures 34 and 35 of [2j. These Figures shows that the spatial variation in effective 
properties correlates well with the variation in volume fraction. With this in mind it appears that 
we should likely choose uj* to be of the size 200 fim. 

We conclude this section by considering a periodic heterogeneous medium of fixed period length 
given by e = \jl > where £ is a fixed positive integer. Here we introduce a method for approxi- 
mation that is derived from the optimal basis associated with the homogenized coefficient obtained 
from periodic homogenization. In what follows we will denote any constant that is independent of 
e and n by C. We write the coefficient describing the periodic medium A e (x) as a rescaling of the 
coefficient of a unit periodic medium, i.e., A e (x) = A(x/e) where A(y) is a coefficient of period one 
for j e II We denote the unit period cell by Q and the homogenized coefficient A is given in 
terms of the periodic corrector matrix P(y), Pij{y) — djw(jj) 1 + Sij where x(u) — (w 1 , u> 2 , w 3 ) is the 
Q periodic solution of 



divA(y)(Vx(y) + I) = 0, for y in 



and 



A a = / P(y)dy. 
Jq 



(6.7) 



(6.8) 



The optimal basis for A e is given in terms of the eigenfunctions ipj of (5.6 1. The optim al basis 
for the homogenized coefficient with A is given in terms of the eigenfunctions ip® of (5.9). We fix 
e = 1/i > and the optimal approximation space is given by the span of the restriction of the 
functions ip\, i = \, ... ,n to uj. In this implementation we introduce the functions <j)\ € Ha' {uj* 



such that 



on duj* , for 



1, 



As before we define the approximation space V e n (w) by 



V?(u) = span{u\ = Pfi, i = 1, 
and state the following approximation theorem. 



,n} 



(6.9) 
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Theorem 6.3. Given any function u € HA'{w*)fR then 



mm 



(6.10) 



where is £/ie n-width associated with A e . Moreover d e n _ 1 is estimated in terms of an easily 

computable quantity 



A e Vui ■ Vui dx 



and the estimate is given by 

\Q2-<r n -i\< ce 1/2 . 

Proof. The theorem is proved by constructing upper bounds on the quantity 



R e = sup <^ inf ||w - w\\ £ j u ) 



(6.11) 



(6.12) 



(6.13) 



From the corrector theory of periodic homogenization it follows from [27 that there exists a constant 
C depending only on ||D 2 (/?° ||l 2 (w*) an d ex < /3 for which 



M-VihHu*) < Ce, 
and since <j)\ £ Ha>(w*)/R with A e G-converging to A H and 

Hence 



(6.14) 

it follows again from [27] that 
(6.15) 

(6.16) 



Now consider it 6 Ha--(lu*)/M. with ||u||£ e ( u «) = 1- For id € V™(uj) there are constants Ci , C2 , . . . , c„ 
such that we can write w — X)"=i c * u i an d we choose these constants c±,C2, . . . ,c n such that 
ip e = ^ji—xCifi gives the optimal approximation to u in the £ e (uj) norm. For this choice one 
has 



X>(tf-* e )k 



< di 



Ce, 



(6.17) 



where the first term on the last line of the ine qualit y follows from the definition of n-width and 
optimal basis and the second term follows from (6.16) and it follows that Rg < d^_ x + Ce. 



We conclude the proof by establishing (6.12). From Theorem 3.1 



K-i) 2 = / A^-Vpldx, 



where we have taken the normalization 



[ A^v t n -Vtf l dx = l. 



(6.18) 



(6.19) 
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On choosing u^ £ V"(uj) we write 



(q:) 2 = [ A e v<-v<dx 

J UJ 



and 



(Q?) 2 - K-i) 2 = / A«(v< + v^)-(v<-vv£)«fc. 

Apriori elliptic estimates show that C = sup e>0 {||u^ + ^£lUe(w*)} < 00 an d 

KQD 2 - K-i) 2 | < cJjf A«(Vu£ - V^) • (Vti{- - Vtf - ) dx 

< c||^-<|U 2(w .) < c e 



(6.20) 



(6.21) 



(6.22) 



where the second to las t inequality follows from Theorem 3.1 and the last inequality follows from 
(6.16). Inequality (6.121 follows noting that 



IQ?-<_il<l(Q?) 2 -(<f„-i)T /2 . 



(6.23) 



A Appendix 



We provide a proof of the Cacciappoli inquality given in Lemma 3.1 We introduce the cut off 
function n £ Cq(uj*) such that < r\ < 1 and r\ = 1 for points inside O and |V7j(x)| < 1/6 for points 
in uj* . Given the function u £ 7i(w*) and since u is A - harmonic we have 



/ AVu ■ V(?7 2 it) dx 

J UJ* 



0. 



(A.l) 



Expanding (A.l) gives 



I {AVu ■ Vu)r? 2 dx = -2 / (nA 1/2 Vu) ■ (uA 1/2 Vn) dx 

J UJ* J UJ* 

<2^ (AVu ■ Vu)r/ 2 dx^j (J (AVr, ■ Vn)u 2 dx 



(A.2) 



\\u\\s(0)<(J (AVu-Vu)rj 2 dx^j < 2 (J (AVrj ■ Vr])u 2 dx\ 
< 272 



1/2 



1/2 2 1/2 

IVTjIVdl) < \\u\\ L 2 {u .y 



(A.3) 



and Lemma |3.1| is proved. 

We now show that the restriction operators introduced in section three are compact. We first 
consider two concentric cubes uj C uj*. The restriction operator P : Ha(ui*)/M. — > Ha(u>)/M. is 
defined by Pu(x) = u(x) for all x £ uj and all u £ Ha(uj*)- 

Lemma A.l. Given any sequence {itn}^Li € Ha(uj*) /K that is bounded in the energy norm over 
(uj*) then one can extract a subsequence that converges in JJ 1 (w) to an element of Ha(uj)/R- 

Proof. We apply the Poincare inequality together wit h the Rellich compactness theorem to extract a 
convergent subsequence in L 2 (uj*). From Lemma 3.1 it now follows that this subsequence is Cauchy 
with respect to the energy norm over uj and the convergence in H (uj) follows. The weak formulation 
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of the boundary value problem together with the strong convergence of the subsequence easily shows 
that the limit function is A-harmonic and the theorem is proved. 

Next we consider two concentric cubes C C w* such that w = Cfl!J and w* n ft have non zero 
volume. Here the side length of C is a and that of ui* is a* = (1 + p)a. The restriction operator 
P : H Ai0 (uj*r}n)/R -> H a . (lo)/R is defined by Pu(x) = u{x) for all x G w and all u G #A,o(w*nfi). 
Here we suppose the boundary of f2 is C . 

Lemma A. 2. Given any sequence G H Ai q(u>* H f2)/K £/ia£ is bounded with respect to the 

energy norm (u>* l~l il) i/ien one can extract a subsequence that converges in H 1 ^) to an element of 
H A , (u)/R. 

Proof. Following section 3 we extend each u n G H\ (lo* n Q)/K as an A-harmonic function across 
dfl onto the set lu* e such that 



l|w n ||.fji(w*) < C||Wn||ffi(a)*nn) 



(A.4) 



where C depends only on <90. Application of Theorem 3.1 gives 



4/31/2 

\S{w) < -^-\\Un\\L2(u* E ) 



(A.5) 



and we deduce that 



4/3^2 

||Wn||fi(w) < C ||wTi||L 2 (Lj*nf2)- 



(A.6) 



With (A.6) in hand we can now proceed as in the proof of Lemma A.l to establish compactness. 
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